% data          :输入数据
% data_len_lim  :数据长度限制,可不填
% Fs            :数据采样率
% spectrum_type :可选项  - "power" : 功率谱计算  -"CUstom" : 频谱计算
% plot_en       :可选项 - "enable": 该函数自动打印图表  -"disable": 该函数不自动打印图表
% plot_tittle   :自动打印图表的标题,仅在plot_en设置为"enable"时有效.

function[spectrum,freq] = plotCustomSpectrum_method1 (data ,data_len_lim, Fs,spectrum_type,plot_en,plot_title)


    %% 将数据数量进行限制,超过这个数据量会导致pspectrum报告FFTW:nomem错误
    if(length(data) < 2400000)
        datalen_max = length(data);
    else
        datalen_max = 2400000;
    end
    if(isempty (data_len_lim))
    else
        datalen_max = data_len_lim;
    end
    timeLimits = [1 datalen_max];
    data_internal = data(timeLimits(1):timeLimits(2));

    %% 计算pspectrum函数参数
    leakage = 0;
    kbeta = 40*(1-leakage);
    enbw = getENBWEstimate(kbeta);
    fresMin = 8*enbw/length(data_internal);
    [spectrum, Fgenerate_wave_ROI] = pspectrum(data_internal,'FrequencyResolution',fresMin,'Leakage',leakage);

    %% 计算输出值
    freq = Fgenerate_wave_ROI * Fs /2 /pi;
    if(spectrum_type == "power")
        spectrum = 10*log10(spectrum);
        y_str = "功率 (dB)";
    else
        spectrum = 20*log10(sqrt(spectrum*2));
        y_str = "幅度 (dB)";
    end
    %% 绘图
    if(plot_en == "enable")
        % ==== 绘图 ====
        figure;
        hold on;

        % ==== 主瓣和噪声划分 ====
        mainlobe_mask = freq <= 20000;  % 主瓣定义为 0~20kHz
        noise_mask    = ~mainlobe_mask;

        % 主瓣用红色线绘制
        plot(freq(mainlobe_mask), spectrum(mainlobe_mask), 'r', 'LineWidth', 1.5);
        % 噪声用蓝色线绘制
        plot(freq(noise_mask), spectrum(noise_mask), 'b--', 'LineWidth', 1);
        % 画辅助线标注主瓣界限
        ylimits = ylim;
        line([22000 22000], ylimits, 'Color', [0.5 0.5 0.5], 'LineStyle', ':');

        % 图形设置
        xlabel('频率 (Hz)');
        ylabel(y_str);
        title(sprintf('%s\n主瓣 vs 噪声对比(主瓣 ≤ 20KHz)', plot_title));
        legend('主瓣频段 (≤20KHz)', '噪声频段 (>20KHz)', 'Location', 'best');
        xlim([0, min(Fs/2, 100e3)]);
        grid on;
    end
end


function enbwValue = getENBWEstimate(kbeta)
%GETENBWESTIMATE Calculate enbw using approximation similar to zoomspectrum engine
%  This function is for internal use only. It may be removed in the future.

%   Copyright 2020 The MathWorks, Inc.
%#codegen

    poly1 = [-7.42100220301407e-05,  0.0010951661775193,  -0.00601630113184729,...
        0.0135330207063678,-0.00553286257981127,   0.00188967508127851,...
        -0.000245892879902476, 1.00000397846243];
    poly2 = [ 5.69809789173636e-06, -0.000159680213659123, 0.00182733327457552,...
        -0.0104618197432362,0.0265970133409903,    0.00682196582785356,...
        -0.0609769138391775,   1.04788712768061];
    poly3 = [ 7.68603379795731e-12, -1.66807848853759e-09, 1.5594811121591e-07,...
        -8.2806698857875e-06,0.000279922339315112, -0.00665206774673593,...
        0.161004938614092,    0.690887010063724];

    if  kbeta <= 2.9
        poly = poly1;
    elseif kbeta <= 4.9
        poly = poly2;
    else
        poly = poly3;
    end
    enbwValue = 0;
    for idx = 1:8
        enbwValue = enbwValue*kbeta + poly(idx);
    end

end